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ABSTRACT 

The development of radiation hydrodynamical methods that are able to follow gas dynamics 
and radiative transfer self-consistently is key to the solution of many problems in numeri- 
cal astrophysics. Such fluid flows are highly complex, rarely allowing even for approximate 
analytical solutions against which numerical codes can be tested. An alternative validation 
procedure is to compare different methods against each other on common problems, in or- 
der to assess the robustness of the results and establish a range of validity for the methods. 
Previously, we presented such a comparison for a set of pure radiative transfer tests (i.e. for 
fixed, non-evolving density fields). This is the second paper of the Cosmological Radiative 
Transfer (RT) Comparison Project, in which we compare 9 independent RT codes directly 
coupled to gasdynamics on 3 relatively simple astrophysical hydrodynamics problems: (5) 
the expansion of an H II region in a uniform medium; (6) an ionization front (I-front) in a 
density profile with a flat core, and (7), the photoevaporation of a uniform dense clump. 
Results show a broad agreement between the different methods and no big failures, indicating 
that the participating codes have reached a certain level of maturity and reliability. However, 
many details still do differ, and virtually every code has showed some shortcomings and has 
disagreed, in one respect or another, with the majority of the results. This underscores the fact 
that no method is universal and all require careful testing of the particular features which are 
most relevant to the specific problem at hand. 

Key words: H II regions — galaxies :high-redshift — intergalactic medium — cosmology: 
theory — radiative transfer — methods: numerical 



1 INTRODUCTION 
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The transfer of ionizing radiation through optically-thick media is a 
key process in many astrophysical phenomena. Some examples in- 
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elude cosmological reionization (e.g. Gnedin 2000; Nakamoto et al. 
2001; Razoumov et al. 2002; Sokasian et al. 2003; Ciardi et al. 
2003; Iliev et al. 2006b; Kohler et al. 2007), star formation (e.g. 
Hosokawa & Inutsuka 2005; Iliev et al. 2006a; Razoumov et al. 
2006; Susa & Umemura 2006; Ahn & Shapiro 2007; Whalen & 
Norman 2008b), radiative feedback in molecular clouds (Mellema 
et al. 2006a; Mac Low et al. 2007; Dale et al. 2007a,b; Krumholz 
et al. 2007; Gritschneder et al. 2009), and planetary nebulae (e.g. 
Mellema et al. 1998; Lim & Mellema 2003). In some of these 
problems, fast, R-type I-fronts predominate. Those fronts propa- 
gate faster than the hydrodynamic response of the gas, so gas mo- 
tions do not affect the I-front evolution. In these cases the radiative 
transfer could be done on a fixed density field (or a succession of 
such fields), and dynamic coupling to the gas is generally not re- 
quired. However, the majority of astrophysical and cosmological 
applications involve slow, D-type I-fronts (or a combination of R- 
type and D-type, as we describe in detail in section 3.1), so the 
radiative transfer and gasdynamics should be directly coupled and 
evolved simultaneously. Until recently, self-consistent radiation hy- 
drodynamical codes for radiative transport have been rare, but this 
unsatisfactory situation is now rapidly changing due to the devel- 
opment of a number of such codes using a variety of numerical 
approaches. 

A number of radiative transfer methods have been developed 
in recent years, both stand-alone and coupled to hydrodynamics. 
High computational costs necessitate the usage of various approx- 
imations. Thus, it is of prime importance to validate the numerical 
methods developed and to evaluate their reliability and accuracy. 
Tests with either exact or good approximate analytical solutions 
should always be the first choice for code testing. Extensive test 
suites of radiation hydrodynamical I-front transport in a variety of 
stratified media with good approximate analytical solutions do exist 
(Franco et al. 1990; Whalen & Norman 2006) and are stringent tests 
of coupling schemes between radiation, gas, and chemistry. How- 
ever, an alternative and complementary approach is to compare a 
variety of methods on a set of well-defined problems in astrophys- 
ical settings. This is the approach we have taken in this project. 

Our aim is to determine the type of problems the codes are 
(un)able to solve, to understand the origin of any differences in- 
evitably found in the results, to stimulate improvements and fur- 
ther developments of the existing codes and, finally, to serve as a 
benchmark for testing future algorithms. All test descriptions, pa- 
rameters, and results can be found at the project website: 
http : // www.cita.utor onto. ca/ ^ iliev /rtwiki/doku.php . 

The first paper of this comparison project discussed the re- 
sults from fixed density field tests (Iliev & et al. 2006, hereafter 
Paper I), i.e. without any gas evolution. We found that all partic- 
ipating codes are able to track I-fronts quite well, within ~10% 
of each other. Some important differences also emerged, especially 
in the derived temperatures and spectral hardening. We found that 
some of these differences were due to variations in microphysics 
(chemical reaction rates, heating/cooling rates and photoionization 
cross-sections), while others were due to the method itself, e.g. how 
the energy equation is solved, how many frequency bins are used 
for the spectral evolution, etc. We concluded that the tested radia- 
tive transfer methods are producing reliable results overall, but that 
not all methods are equally appropriate for any given problem, es- 
pecially in cases when obtaining precise temperatures and spectral 
features is important. 

We now extend our previous work by considering a set of radi- 
ation hydrodynamical tests. In the spirit of Paper I, we have chosen 
a set of test problems which are relatively simple, so as to be most 
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Figure 1. Legend for the line plots. 

inclusive given the current limitations of the available codes (e.g. 
1-D or 2-D vs. 3-D codes). At the same time, our tests consider 
problems of astrophysical importance, and cover a wide variety of 
situations that test the attributes of each method, including its ra- 
diative and hydrodynamic components and their coupling. 

The efficiency, optimization and performance of the codes 
are very important, especially for the most complex and 
computationally-intensive problems. However, there are a number 
of complications, which we discussed in Paper I, preventing us 
from doing such testing in a meaningful way at present. We there- 
fore leave it for future work. 

All test results for this study had to be supplied on a regular 
Cartesian grid of 128^ computational cells. This relatively modest 
resolution was chosen in the interests of inclusivity, so that even 
codes which are not yet fully optimized in terms of either compu- 
tations or memory can participate in the comparison. We note that 
production runs at present are typically run at 256^ or better resolu- 
tion. Codes which utilize Adaptive Mesh Refinement (AMR) grids 
or particles have been requested to run the problem at the resolu- 
tion which approximates as closely as possible the fixed-grid one 
for fair comparison. Their results have then been interpolated onto 
a regular grid for submission. 



2 THE CODES 

In this section we briefly describe the nine radiative transfer codes 
participating in this stage of the comparison project, with refer- 
ences to more detailed method papers if available. Details of the 
codes and their basic features and methods are summarized in Ta- 
ble 1. Figure 1 provides a legend allowing the reader to identify 
which line corresponds to which code in the figures throughout the 
paper. The images we present are identified in the corresponding 
figure caption. 

2.1 Capreole+C^-Ray and TVD+C^-Ray (G. Mellema, I. 
Iliev, P. Shapiro, M. Alvarez) 

C^-Ray (Mellema et al. 2006b) is a grid-based short characteristics 
(e.g. Raga et al. 1999) ray-tracing code which is photon-conserving 
and causally traces the rays away from the ionizing sources up to 



Table 1. Participating codes and their current features. 
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Code 



Grid 



Parallelization 



hydro method 



rad. transfer method 



Capreole+C^ -Ray 

TVD+C^-Ray 

HART 

RSPH 

ZEUS-MP 

RHID 

Coral 

LICORICE 

Flash-HC 

Enzo-RT 



fixed shared/distributed 

fixed shared/distributed 

AMR shared/distributed 

particle-based distributed 

fixed distributed 

sph. Lagrangian no 

AMR no 

AMR shared 

AMR distributed 

fixed distributed 



Eulerian, Riemann solver 

Eulerian, TVD solver 

Eulerian, Riemann solver 

SPH 

Eulerian 

Lagrangian 

Eulerian, flux- vector splitting 
SPH 

Eulerian, PPM 
Eulerian, PPM 



short-characteristics ray-tracing 
short-characteristics ray-tracing 
Eddington tensor moment 
long-characteristics ray-tracing 
3-D ray-tracing 
1-D ray-tracing 

short-characteristics ray tracing 
Monte-Carlo ray-tracing 
Hybrid characteristics ray-tracing 
Flux-limited diffusion 



each cell. Explicit photon-conservation is assured by taking a finite- 
volume approach when calculating the photoionization rates, and 
by using time- averaged optical depths. The latter property allows 
for integration time steps that are much larger than the ionization 
time scale, which results in a considerable speed-up of the calcula- 
tion and facilitates the coupling of the code to gasdynamic evolu- 
tion. The code and the various tests performed during its develop- 
ment are described in detail in Mellema et al. (2006b). 

The frequency dependence of the photoionization rates and 
photoionization heating rates are dealt with by using frequency- 
integrated rates, stored as functions of the optical depth at the ion- 
ization threshold. In its current version the code includes only hy- 
drogen and no helium, although it could be added in a relatively 
straightforward way. 

The transfer calculation is done using short characteristics, 
where the optical depth is calculated by interpolating values of 
grid cells lying along the line-of-sight towards the source. Because 
of the causal nature of the ray-tracing, the calculation cannot eas- 
ily be parallelized through domain decomposition. However, us- 
ing OpenMP and MPI the code is efficiently parallelized over the 
sources and grid octants (Iliev et al. 2008b). The code has been 
applied for large-scale simulations of cosmic reionization and its 
observability (Iliev et al. 2006b; Mellema et al. 2006c; Iliev et al. 
2007a,b; Holder et al. 2007; Dore et al. 2007; Iliev et al. 2008a) on 
grid sizes up to 406^ and up to ~ 10^ ionizing sources, running on 
up to 10,240 computing cores. 

There are ID, 2D and 3D versions of the code that are avail- 
able. It was developed to be directly coupled with hydrodynamics 
calculations. The large time steps allowed for the radiative trans- 
fer enable the use of the hydrodynamic time step for evolving the 
combined system. The C^-Ray radiative transfer and nonequilib- 
rium chemistry code has been coupled to several different gasdy- 
namics codes, utilizing both fixed and adaptive grids. The tests in 
this project were mostly performed with the version coupled to the 
hydrodynamics code Capreole developed by Garrelt Mellema and 
based on Roe's approximate Riemann solver. The first gasdynamic 
application of our code is presented in Mellema et al. (2006a). Ad- 
ditionally, one of the tests has also been run with C^-Ray coupled 
to a different hydro solver, namely the TVD method of Trac & Pen 
(2004) (see Test 6 below). 



2.2 Hydrodynamic Adaptive Refinement Tree (HART) (N. 
Gnedin, A. Kravtsov) 

The Hydrodynamic Adaptive Refinement Tree (HART) code is an 
implementation of the AMR technique and uses a combination of 
multi-level particle-mesh and shock-capturing Eulerian methods 
for simulating the evolution of the dark matter particles and gas, 
respectively. High dynamic range is achieved by applying adaptive 
mesh refinement to both gas dynamics and gravity calculations. 

The code performs refinements locally on individual cells, and 
cells are organized in refinement trees (Khokhlov 1998). The data 
structure is designed both to reduce the memory overhead for main- 
taining a tree and to fully eliminate the neighbor search required for 
finite-difference operations. All operations, including tree modifi- 
cations and adaptive mesh refinement, can be performed in parallel. 
The advantage of the tree-based AMR is its ability to control the 
computational mesh on the level of individual cells. This results in 
a very efficient and flexible (and thus highly adaptive) refinement 
mesh which can be easily built and modified and, therefore, effec- 
tively match the complex geometry of cosmologically interesting 
regions: filaments, sheets, and clumps. Several refinement criteria 
can be combined with different weights allowing for a flexible re- 
finement strategy that can be tuned to the needs of each particular 
simulation. The adaptive refinement in space is accompanied by a 
temporal refinement (smaller time steps on meshes of higher reso- 
lutions). 

The ART code was initially developed by A. Kravtsov in col- 
laboration with A. A. Klypin and A. M. Khokhlov (Kravtsov et al. 
1997; Kravtsov 1999; Kravtsov et al. 2002). N. Gnedin joined 
the HART code development team in the spring of 2003 and has 
adopted the OTVET algorithm for modeling 3D radiative transfer 
for the ART mesh structure and implemented a non-equilibrium 
chemical network and cooling (e.g. Gnedin et al. 2008). 



2.3 RSPH (H. Susa, M. Umemura, D. Sato) 

The Radiation-SPH (RSPH) scheme is specifically designed to in- 
vestigate the formation and evolution of first-generation objects at 
z ^ 10, where the radiative feedback from various sources plays 
important roles. The code can compute the fraction of chemical 
species e~, H+, H, H~, H2, and H^ by fully implicit time integra- 
tion. It also can deal with multiple sources of ionizing radiation, as 
well as with Lyman- Werner band photons. 

Hydrodynamics is calculated by the smoothed particle hydro- 
dynamics (SPH) method. It uses the version of SPH by Umemura 
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(1993) with the modification according to Steinmetz & Mueller 
(1993), and adopts the particle resizing formalism by Thacker et al. 
(2000). The present version does not use the so-called entropy for- 
malism (Springel & Hernquist 2002). The non-equilibrium chem- 
istry and radiative cooling for primordial gas are calculated using 
the code developed by Susa & Kitayama (2000), where H2 cooling 
and reaction rates are taken from Galli & Palla (1998a). 

As for the photoionization process, the on-the-spot approxi- 
mation is employed (Spitzer 1978), meaning that the transfer of 
ionizing photons directly from the source is solved, but diffuse pho- 
tons are not transported. Instead, it is assumed that recombination 
photons are absorbed in the same zone from which they are emit- 
ted. Due to the absence of the source term in this approximation, 
the radiation transfer equation becomes very simple. Solving the 
transfer equation reduces to the easier problem of assessing the op- 
tical depth from the source to every SPH particle. 

The optical depth is integrated utilizing the neighbour lists 
of SPH particles. It is similar to the code described in Susa 
& Umemura (2004), but can now also deal with multiple point 
sources. In the new scheme fewer grid points are created on the 
light ray than in its predecessor. Instead, just one grid point per SPH 
particle is created in the particle's neighborhood. The 'upstream' 
particle for each SPH particle on its line of sight to the source is 
then found. Then the optical depth from the source to the SPH par- 
ticle is obtained by summing up the optical depth at the 'upstream' 
particle and the differential optical depth between the two particles. 

The code is parallelized with the MPI library. The computa- 
tional domain is divided by the Orthogonal Recursive Bisection 
method. The parallelization method for radiation transfer is sim- 
ilar to the Multiple Wave Front method developed by Nakamoto 
et al. (2001) and Heinemann et al. (2006), but it is adapted to fit the 
SPH code as described in (Susa 2006). 

The code computes self-gravity using a Barnes-Hut tree, 
which is parallelized as well. A Tree-GRAPE version of the code 
has also been developed. This code has been applied to radiative 
feedback in primordial star formation (Susa & Umemura 2006; 
Susa 2007; Hasegawa et al. 2009), as well as the regulation of 
star formation in forming galaxies by ultraviolet background (Susa 
2008). 



2.4 ZEUS-MP (D. Whalen, J. Smidt, M. Norman) 

ZEUS -MP solves explicit finite-difference approximations to Eu- 
ler's equations of fluid dynamics self-consistently with a 9-species 
primordial gas reaction network (H, H^, He, He^, He^^, H~, H2, 
H J and e) and ray-tracing radiative transfer, which is used to com- 
pute the radiative rate coefficients required by the network and the 
gas energy equation. Our method is described in detail elsewhere 
(Whalen & Norman 2006, 2008b); here, we review multifrequency 
upgrades to the radiative transfer and improvements to the subcy- 
cling scheme (Whalen & Norman 2008a). 

The ZEUS-MP RT module evaluates radiative rate coefficients 
by solving the static equation of transfer in flux form. To obtain the 
total rate coefficient k for a zone we sum the ku computed for a 
given binned photon emission rate over all energies by looping the 
solution to the transfer equation over them. In tests spanning 40 
to 2000 energy bins, good convergence is found with 120 bins, 40 
bins spaced evenly in energy from 0.755 eV to 13.6 eV and 80 bins 
that are logarithmically- spaced from 13.6 eV to 90 eV. 

Successive updates to the reaction network and gas energy are 
performed over the minimum of the chemical time 



_ ne+O.OOlng 

tchem — U.i : . (l; 

Tie 

and the photoheating/cooling time 

t^, = 0.1^^^ (2) 

^ht/ cool 

until the larger of these two times has been crossed, at which point 
full hydrodynamical updates of gas densities, energies, and veloc- 
ities are performed. These times are global minima for the entire 
grid. Chemical times are defined in terms of electron flow to ac- 
commodate all chemical processes rather than just ionizations or 
recombinations. Adopting the minimum of the two times for chem- 
istry and gas energy updates enforces accuracy in the reaction net- 
work when tchem becomes greater than the (in relic H II regions, 
for example). 

ZEUS-MP is now fully parallelized for three-dimensional ap- 
plications. We have updated the H and He recombination and cool- 
ing rates responsible for some minor departures between ZEUS- 
MP and the other codes in Paper I in the temperature structure 
of H II regions, and now use the most recent data from Hummer 
(1994) and Hummer & Storey (1998). Our code has been validated 
with stringent tests of R-type and D-type I-fronts in a variety of 
stratified media (Franco et al. 1990; Whalen & Norman 2006) and 
applied to both cosmological and astrophysical problems, such as 
the breakout of UV radiation from primordial star-forming clouds 
(Whalen et al. 2004), the formation of dynamical instabilities in 
galactic H II regions (Whalen & Norman 2008b), the circumstel- 
lar environments of gamma-ray bursts (Whalen et al. 2008b), the 
photoevaporation of cosmological minihalos by nearby primordial 
stars (Whalen et al. 2008a), and Pop III supernovae explosions in 
cosmological H II regions (Whalen et al. 2008c). 

2.5 RHID (K. Ahn, P. Shapiro) 

RHID is a ID, Lagrangian, spherically- symmetric, radiation- 
hydrodynamics code for a two-component gas of baryons and col- 
lisionless dark matter coupled by gravity (Ahn & Shapiro 2007). 
For the baryonic component, the Euler equations and the equation 
of state are solved, together with multi-frequency, multi-species 
radiative transfer equations and a reaction network with nine pri- 
mordial species (H, H+, He, He+, He++, H", H2, H+ and e). 
Dark matter dynamics, governed by the coUisionless Boltzmann 
equations, takes a simplified form in spherical symmetry. The code 
solves an effective set of Euler equations for a dark matter fluid, 
based upon the "fluid approximation" of dark matter dynamics for 
a spherically symmetric system with an isotropic velocity disper- 
sion, derived and justified elsewhere (Ahn & Shapiro 2005). These 
effective Euler equations are identical to those for an inviscid, ideal 
gas with a ratio of specific heats 7 = 5/3. 

The Euler equations are solved using the so-called "leap-frog" 
method, where the Lagrangian position (radius) and velocity (ra- 
dial velocity) are staggered in time to achieve a second-order accu- 
racy in time steps, both for baryonic and dark matter fluid. The 
usual artificial viscosity scheme is used to capture shocks. We 
typically adopt a few thousand uniformly spaced bins in radius. 
Non-equilibrium rate equations for the nine primordial species are 
solved using the backward differencing scheme of Anninos et al. 
(1997). For H~ and HJ, due to their relatively fast reaction rates, 
the equilibrium values may be used. 

Radiative transfer is performed by ray-tracing, taking account 
of the optical depth to bound-free opacity of H I, He I, He II, H~, 
and H2, as well as bound-free and dissociation opacity of H^. The 
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optical depth to the Lyman- Werner band photons of H2, which are 
capable of dissociating H2, is treated using a pre-calculated self- 
shielding function by Draine & Bertoldi (1996), which is deter- 
mined by the H2 column density and gas temperature. Diffuse flux 
is not explicitly calculated, but is accounted for implicitly by adopt- 
ing case B recombination rates. The radiative reaction rates are 
calculated using a photon-conserving scheme, which enables the 
code to treat optically- thick shells (e.g. Razoumov & Scott 1999; 
Abel et al. 1999a). A wide range of radiation frequency (energy), 
hiy ^ [0.7 — 7000] eV, is covered by a few hundred, logarithmi- 
cally spaced bins, together with additive, linearly spaced bins where 
radiative cross sections change rapidly as frequency changes. For 
each frequency and species, the corresponding radiative reaction 
rate is calculated, then summed over frequency to obtain the net 
radiative reaction rate. 

The radiative transfer scheme is able to treat 1) an internal 
point source, 2) an external, radially-directed source, and 3) an ex- 
ternal, isotropic background. The transfer for (1) and (2) is ID, per- 
formed along the radial direction only. For (3), the transfer is 2D in 
nature, and at each point the mean intensity is required to calculate 
the radiative rates, which involves an angle integration. The radia- 
tive transfer calculation is performed for each pre-selected angle (6, 
measured from the radial direction), and then the angle integral is 
calculated using the Gaussian quadrature method. 

The code adopts a very stringent time step criterion for accu- 
racy. The minimum of dynamical, sound-crossing, cooling/heating, 
and species change time scales, which is multiplied by a coefficient 
smaller than unity 0.1), is chosen as the time step. All the Euler 
equations and rate equations are solved with this time step, which 
makes the whole calculation self-consistent. This code has been 
tested extensively and used to study the radiative feedback effects 
by the first stars on their nearby minihalos (Ahn & Shapiro 2007). 

2.6 Coral (I. Iliev, A. Raga, G. Mellema, P. Shapiro) 

CORAL is a 2-D, axisymmetric Eulerian fluid dynamics AMR 
code (see Mellema et al. 1998; Shapiro et al. 2004, and refer- 
ences therein for detailed description). It solves the Euler equations 
in their conservative finite-volume form using the second-order 
method of van Leer flux- splitting, which allows for correct and pre- 
cise treatment of shocks. The grid refinement and de-refinement 
criteria are based on the gradients of all code variables. When the 
gradient of any variable is larger than a pre-defined value the cell is 
refined, and when the criterion for refinement is not met the cell is 
de-refined. 

The code follows, by a semi-implicit method, the non- 
equilibrium chemistry of multiple species (H, He, C II- VI, N I- 
VI, O I- VI, Ne I- VI, and S II- VI) and the corresponding cooling 
(Raga et al. 1997; Mellema et al. 1998), as well as Compton cool- 
ing. The photoheating rate is the sum of the photoionization heat- 
ing rates for H I, He I and He II. For computational efficiency all 
heating and cooling rates are pre-computed and stored in tables. 
The microphysical processes - chemical reactions, radiative pro- 
cesses, transfer of radiation, heating and cooling - are implemented 
though the standard approach of operator- splitting (i.e. solved at 
each time- step, side-by-side with the hydrodynamics and coupled 
to it through the energy equation). The latest versions of the code 
also include the effects of an external gravity force. 

Currently the code uses a black-body or power-law ionizing 
source spectrum, although any other spectrum can be accommo- 
dated. Radiative transfer of the ionizing photons is treated explicitly 
by taking into account the bound-free opacity of H and He in the 



photoionization and photoheating rates. The photoionization and 
photoheating rates of H I, He I and He II are pre-computed for the 
given spectrum and stored in tables vs. the optical depths at the ion- 
izing thresholds of these species, which are then used to obtain the 
total optical depths. The code correctly tracks both fast (by evolv- 
ing on an ionization timestep. At ~ riu/nu) and slow I-fronts. 

The code has been tested extensively and has been applied to 
many astrophysical problems, e.g. photoevaporation of clumps in 
planetary nebulae (Mellema et al. 1998), cosmological minihalo 
photoevaporation during reionization (Shapiro et al. 2004; Iliev 
et al. 2005), and studies of the radiative feedback from propagat- 
ing ionization fronts on dense clumps in damped Lyman-a systems 
(Iliev et al. 2006a). 

2.7 LICORICE: Line COntinuum Radiative tranfer 

Integrated Computing Engine (S. Baek, B. Semelin, F. 

Combes) 

The LICORICE code has three main components: TreeSPH to 
compute gravity and hydrodynamics, continuum radiative trans- 
fer with hydrogen and helium ionization physics, and Lyman-alpha 
line transfer. The latter is not relevant to this comparison and has 
been described elsewhere. The ionizing continuum transfer has 
been described in details in Baek et al. (2009). 

The current version of LICORICE does not include H2 for- 
mation, or diffuse radiation from recombinations, but they will be 
incorporated in the future. LICORICE uses SPH particles for the 
gas dynamics and an adaptive grid for the radiative transfer. Physi- 
cal quantities are interpolated from one to the other as required. 

The fluid dynamics are followed using a TreeSPH method. 
The implementation is described in detail in Semelin & Combes 
(2002) and Semelin & Combes (2005). Since there are many va- 
rieties of SPH, we summarize the main features of our algorithm 
here. We use a spherically-symmetric spline-smoothing kernel and 
50 neighbours to compute the SPH quantities using an arithmetic 
average between the neighbours of the smoothing length h and the 
simple viscosity scheme by Monaghan (1992). 

For the tests in this paper we implemented transmissive 
boundary conditions. This was achieved as follows: for each SPH 
particle within a distance of the simulation box boundary smaller 
than its smoothing length /i, we create a symmetrical 'ghost' par- 
ticle on the other side of the boundary. All physical quantities for 
this ghost particle are equal to those of the initial particle, including 
the velocity. The ghost particles are used as neighbours to compute 
the SPH quantities of real particles. The ghost particles are erased 
and recreated at each time step. 

The continuum radiative transfer is solved using a Monte 
Carlo approach similar to the one employed in the CRASH code 
(Maselli et al. 2003). Here we summarize only the differences be- 
tween LICORICE and CRASH. We compute the gas density at 
each particle's position with the SPH smoothing kernel, and phys- 
ical quantities such as ionization fraction and temperature are up- 
dated according to these particle densities. The density field is gen- 
erally smooth, but may sometimes show spurious fluctuations if the 
particle number density changes sharply. This is a well known but 
unavoidable problem with SPH. 

The radiation field is discretized into photon packets and prop- 
agated through cells along directions chosen at random. The cells 
form an adaptive grid which is derived from the tree structure of 
the particle distribution. Our adaptive grid is built to keep the num- 
ber of particles in each cell within a given range (1 to 8 and 1 to 1 
ranges have been used). This yields greater resolution in the denser 
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regions. The adaptive grid also requires fewer cells than a fixed 
grid to best sample a given inhomogeneous particle distribution, 
thus saving both memory and CPU time. 

The time step for updating physical quantities within a cell 
is also adaptive. We update the physical quantities for all cells and 
particles after the propagation of the number of photon packets cor- 
responding to an integration time dt. However, if the number of ac- 
cumulated photons in a cell during this integration time is greater 
than a pre-set limit {e.g. 10% of the total number of neutral hydro- 
gen atoms in the cell), we update the physical quantities in this cell 
with a time step dt' < dt corresponding to the time elapsed since 
the last update. 

The test results are interpolated from the particle distribu- 
tion onto the 128^ uniform Cartesian grids required in this study. 
Currently, the dynamical part of the code is parallelized for both 
shared and distributed memory architectures using OpenMP and 
MPI, while the radiative transfer is parallelized with OpenMP only. 
The code can now handle 256^ particles, to be increased to 512^ 
in the near future. We note that compared to a uniform grid with 
the same number of cells, the SPH Lagrangian approach results in 
higher resolution in the dense regions, but lower resolution in more 
diffuse regions. 

2.8 Flash-HC: Hybrid Characteristics (T. Theuns, M. 
Raicevic, E.-J. Rijkhorst) 

The Hybrid Characteristics (HC) method (Rijkhorst 2005; Ri- 
jkhorst et al. 2006) is a three-dimensional ray-tracing scheme for 
parallel AMR codes. It combines elements of long and short char- 
acteristics, using the precision and parallelizability of the former 
with efficient execution through interpolation of the latter. It has 
been implemented into the Flash-HC AMR code (Fryxell et al. 
2000), enabling simulations of radiation hydrodynamics problems 
with point sources of radiation. The public version of the Flash 
code (which does not currently include this radiative transfer mod- 
ule) can be downloaded from 

http : // flash.uchicago.edu/website/ home/ . 

The block-structured AMR grid used in Flash-HC is dis- 
tributed over processors using a space-filling curve. Parallel ray 
tracing requires each ray to be split in the independent sections 
where the ray traverses the blocks held by a given processor. First, 
every processor traces rays on its local blocks in directions which 
start from the source, and end in the corners of each cell on the faces 
of the (cubic) block. Since rays cross several blocks, interpolation 
is used to assemble a ray from local block contributions. However, 
because some of these blocks will be held by other processors, local 
column densities need to be exchanged in one global communica- 
tion. Note that only face values are exchanged. Finally, local and 
imported column densities are combined using interpolation to as- 
semble the complete ray. At the end of this parallel operation, each 
cell has the total column density to the source along a ray that tra- 
verses all intervening cells at the full resolution of the AMR grid. 
Interpolation coefficients are chosen such that the exact solution for 
the column density is obtained for a uniform density distribution. 
Even in a non-uniform density distribution, for example the 
differences between the value of the correct column density and 
that obtained using HC is typically less than half a percent. 

Recent improvements introduced since Paper I include the im- 
plementation of a fully photon conserving chemistry solver, taking 
into account the effects of both spatial and temporal discretization 
(Abel et al. 1999b; Mellema et al. 2006b). This implementation 
employs the Livermore Solver for Ordinary Differential Equations 



(LSODE Hindmarsh 1980), which, although more computation- 
ally intensive than the original solver used in DORIC (Frank & 
Mellema 1994), eliminates the need for an independent radiative 
transfer time step irrespective of the ionization front type and guar- 
antees correct front positions and ionization heating. Details of the 
scheme will be presented elsewhere. Additional functionality al- 
lows for a radiation source outside of the computational volume, a 
feature used in Test 7 to approximate a parallel ionization front. 

The parallel scaling of HC was examined in Rijkhorst (2005) 
and Paper I; the algorithm scales well for 100 processors on a 
SGI Altix, and ~ 1000 processors on a IBM Blue Gene/L sys- 
tems. The algorithm scales linearly with the number of sources. 
The photon-conserving RT and chemistry upgrades should not af- 
fect HC's scaling. 

2.9 Enzo-RT (D.R. Reynolds, M.L. Norman, J.C. Hayes, P. 
Paschos) 

The Enzo-RT code is an extension to the freely-available Enzo 
code^ that self-consistently incorporates coupled radiation trans- 
port and chemical ionization kinetics within Enzo's formulation 
for cosmological hydrodynamics on AMR meshes (Norman et al. 
2007). In Enzo-RT we approximate the radiation transport pro- 
cesses using a single integrated radiation energy density in each 
spatial cell that is propagated with flux-limited diffusion on a finite- 
volume mesh. The radiation field is implicitly coupled in time 
to a multi- species chemical reaction network. This implicit radi- 
ation chemistry system is then coupled in an operator- split fashion 
with Enzo's cosmological hydrodynamics solver, which utilizes the 
Piecewise Parabolic Method for the advection of matter and gas 
energy (Colella & Woodward 1984). The coupled algorithm, along 
with a suite of verification tests, is fully described in Reynolds et al. 
(2009). 

The frequency dependence of the photoionization rates is 
treated by integrating a prescribed radiation frequency spectrum, 
typically chosen to be either monochromatic, blackbody, or a 
(^) ^ power law. This integration is performed upon initializa- 
tion of the solver and the integrated rates are re-used throughout 
the simulation. In the current version of the code, only a single ra- 
diation profile is allowed, although this formulation may be easily 
extended to allow for multifrequency calculations. 

The solver for propagating radiation throughout the domain 
follows a standard flux-limited diffusion model, in which the radi- 
ation flux F is approximated by 

F^-i^V^;. (3) 

a 

Here E is the radiation energy density, and the flux-limiter D 
smoothly connects the limiting cases of (nearly) isotropic and free- 
streaming radiation: 



D{E) = diag{Di{E),D2{E),D3{E)) 
where 

c{2kt + Ri) 



Di{E) = 



z = l,...,3, 



(4) 



(5) 



and Ri = \diE\/E, c is the speed of light, and k,t is the opacity. 

The coupled implicit radiation-chemistry system further in- 
cludes a gas energy feedback field, which allows us to self- 
consistently heat and cool the gas in an operator-split manner, 

^ http://lca.ucsd.edu/portal/software/enzo 
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Figure 2. Test 5 (H II region expansion in an initially-uniform gas): Images of the H I fraction, cut through the simulation volume at coordinate z = at time 
t = 100 Myr for (left to right and top to bottom) Capreole+C^ -Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



capturing all of the stiff components involved in radiation trans- 
port, primordial chemistry and thermal heating/cooling in a tightly- 
coupled implicit system. Enzo's explicit Eulerian hydrodynamics 
solver and its parallel implementation have been exhaustively de- 
scribed elsewhere (Norman et al. 2007). Parallelism of the cou- 
pled implicit system follows a standard domain-decomposition 
approach and is solved using state-of-the-art Newton- Krylov- 
Multigrid solvers (Knoll & Keyes 2004), potentially allowing scal- 
ability of the algorithm to up to tens of thousands of processors. 

While Enzo allows for spatial adaptivity through structured 
adaptive mesh refinement (SAMR), our initial implementation of 
Enzo-RT is currently limited to uniform grids in 1-, 2- or 3- 
dimensions, although their upgrade to AMR is under development. 
Extensions of this approach to variable Eddington tensors, multi- 
group flux-limited diffusion, or multigroup variable Eddington ten- 
sors are easily accommodated within our implicit formulation and 
are planned as future extensions. One benefit is that the timestep 
is independent of grid resolution, at least for the radiation solve. 
Another advantage of our approach is that by defining radiation 
as a field variable, scalability with respect to the number of point 
sources ceases to be an issue. Instead, scalability is dictated by the 
underlying linear system solver, which for the case of multigrid is 
optimal. 



3 RADIATION HYDRODYNAMICS TESTS: 
DESCRIPTION 

For simplicity and inclusivity (since currently not all codes have 
implemented helium or metals chemistry and cooling) all tests as- 
sume the gas to be composed of pure hydrogen. 



3.1 Test 5: Classical H II Region Expansion 

Test 5 is the classical problem of the expansion of an I-front 
due to a point source in an initially uniform-density medium. In 
general, I-fronts are classified according to their speed with re- 
spect to the gas and the change in gas density through the I-front 
(c.f. Kahn & Dyson 1965; Spitzer 1978). There are two critical 
speeds: R-critical, defined sls vr = 2cs,i,2, and D-critical, given 
by VD = Cs,i,2 - {cij^2 - ~ cj^j^i/(2cs,i,2), where 

Cs,i,i — (pi/pi)"^^^ and Cs, 1,2 — (^2/^2)"^^^ are the isothermal 
sound speeds in the gas ahead of and behind the I-front, respec- 
tively. Note that in the test the gas is not assumed to be isothermal. 
The velocity of the I-front is given by the jump condition vi — F/n 
(which guarantees photon conservation), where n is the number 
density of the neutral gas entering the front and F is the flux of 
ionizing photons at the I-front transition (which is attenuated due 
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to absorptions in the gas on the source side). We note that this jump 
condition is modified significantly for I-fronts moving with rela- 
tivistic speeds with respect to the gas (Shapiro et al. 2006). This can 
occur in a number of astrophysical and cosmological environments. 
However, we do not consider such cases here since currently few 
radiative transfer codes (and no radiation hydrodynamics codes, to 
our knowledge) are able to handle such relativistic I-fronts. 

When VI vr (e.g. close to the source, where the flux F is 
large) the I-front is R-type (R-critical when vi = vr). R-type I- 
fronts always move supersonically with respect to the neutral gas 
ahead, while with respect to the ionized gas the front can move ei- 
ther subsonically (strong R-type, highly compressive, but generally 
irrelevant to H II regions since it means that the isothermal sound 
speed behind the front is lower than the one ahead of it), or super- 
sonically (weak R-type, resulting in only slight compression of the 
gas moving through the front). When vi ^ vd, the I-front is D- 
type (D-critical in the case that vi = vd)- The gas passing through 
this type of I-front always expands, and the front is subsonic with 
respect to the gas beyond. With respect to the ionized gas, the I- 
front can again be either supersonic (strong D-type), or subsonic 
(weak D-type). When vd < vi < vr (sometimes referred to as 
an M-type I-front) the I-front is necessarily led by a shock which 



compresses the gas entering the I-front sufficiently to slow it down 
and guarantee that vi ^ vd. 

In a static medium with number density nu and constant ion- 
ized gas temperature T, the evolution of the I-front radius n and 
velocity vi for a point source emitting Nj ionizing photons per 
second are given by 

1/3 



n 



vi ^ 



where 







n 
rs 



- exp(-t/trec) 

exp (— t/tr 



3tr, 



[1 - exp(-t/trec)] 



2/3 



47ras(T)n^ 



1/3 



(6) 
(7) 

(8) 



the Stromgren radius (assuming full ionization), which is reached 
when the number of recombinations in the ionized volume per unit 
time exactly balances the number of ionizing photons emitted by 
the source per unit time. This final static stage is commonly referred 
to as Stromgren sphere. The recombination time is given by 



[aB{T)nu] 



(9) 



Here a b (T) is the case B recombination coefficient of hydrogen at 
temperature T in the ionized region. 
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Figure 4. Test 5 (H II region expansion in an initially-uniform gas): Images of the temperature, cut through the simulation volume at coordinate z ■ 
t = 100 Myr for (left to right and top to bottom) Capreole+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



: at time 



In reality, the ionized gas is not static and its much higher pres- 
sure than that of the ambient medium causes it to expand outward 
beyond the Stromgren radius. Analytical models predict that in this 
phase the I-front radius evolves approximately according to (c.f. 
Spitzer 1978) 



ri 







1 + 



7cst 
4^ 



4/7 



(10) 



where rs is the Stromgren radius and Cs is the sound speed in the 
ionized gas. The expansion finally stalls when a pressure equilib- 
rium is reached. The predicted final H II region radius is 



TeJ 



rs, 



(11) 



where T is the temperature inside the H II region and Te is the 
external temperature. In reality, the evolution is more complicated, 
with non-uniform temperatures inside the H II region, broadened 
I-fronts due to pre-heating by energetic photons, etc. Furthermore, 
equation 10 describes correctly only in the purely pressure-driven, 
late-time evolution, but not the transition from fast, R-type to D- 
type I-front. These analytical solutions should therefore only be 
considered to be guidelines for the expected behaviour, not as exact 
solutions for this problem. 



The numerical parameters for Test 5 are as follows: computa- 
tional box size L = 15 kpc, initial gas number density uh = 10~^ 
cm~^, initial ionization fraction x = 0, constant ionizing pho- 
ton emission rate iV^ = 5 x 10^^ s~^, initial gas velocity zero 
and initial gas temperature Te = 100 K. The radiation source is 
at the (xs^Vs^Zs) = (0,0,0) corner of the computational box. 
For reference, if we assume that the temperature of the ionized 
gas is T = 10^ K, and that the recombination rate is given by 
asiT) = 2.59 x 10"^^ cm^s"\ we find tree = 3.86 x 10^^ s = 
122.4 Myr, rs — 5.4 kpc, and r/ ~ 185 kpc. This rough fi- 
nal pressure-equilibrium radius is thus well outside of our com- 
putational volume, which was instead chosen to resolve the more 
physically-interesting transition from R-type to D-type, which oc- 
curs around r%. Boundary conditions are reflective for the bound- 
aries which contain the origin (where the ionizing source is po- 
sitioned) and transmissive for the other boundaries. The ionizing 
spectrum is that of a 10^ K black body, as expected for a mas- 
sive, metal-free Pop III star. Hydrogen line cooling, recombina- 
tional cooling, and bremsstrahlung cooling are all included, but not 
Compton cooling. The simulation running time is tsim = 500 Myr 
^ 4 tree. Thc rcquircd outputs are the neutral fraction of hydrogen, 
gas pressure, temperature and Mach number on the entire grid at 
t = 10, 30, 100, 200, and 500 Myr, and the I-front position (de- 
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Figure 5. Test 5 (H II region expansion in an initially-uniform gas): Images of the H I fraction, cut through the simulation volume at coordinate z = at time 
t = 500 Myr for (left to right and top to bottom) Capreole+C^ -Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



fined as the position where the neutral fraction is 50%) and I-front 
velocity vs. time along the x-axis. 



3.2 Test 6: H II region expansion in density profile 

Test 6 is the propagation of an I-front created by a point source 
at the center of a spherically-symmetric, steeply-decreasing power- 
law density profile with a small flat central core of gas number den- 
sity no and radius ro : 



nH(r) 



no 



noir/roY 



if r ^ ro 
if r ^ ro 



For a static-density medium the evolution of the I-front within 
the flat-density core is described by equations 6 and 7. In this 
case, if the Stromgren radius associated with the core density no, 
rs,o = [3iV^/(47ras(T)no)]^/^, is smaller than ro, the front will 
come to a halt within the core. If, instead, rs,o > ro, the front 
escapes the core and propagates into the stratified envelope. There- 
after, the I-front position and velocity as a function of time have 
complex analytical forms for an arbitrary source fluxes and densi- 
ties (Mellema et al. 2006b). A simple solution exists for the spe- 
cial case of the central ionizing source rate of photon emission 



= 167rronoas/3, in which case the I-front radius upon leav- 
ing the core is 



n = ro(l + 2t/trec,core)^^^, 



(12) 



where tree, core IS thc recombination time in the core (Mellema et al. 
2006b). Similar solutions exist also when the I-front is moving rel- 
ativistically (Shapiro et al. 2006). 

The propagation of an I-front in r~^ density profiles with full 
gas dynamics does not have an exact analytical solution, but has 
been well studied with both semianalytical and numerical methods 
(Franco et al. 1990). If rs,o < ro then the I-front converts to D- 
type within the core, but starts to re-accelerate upon entering the 
steep density gradient. Numerical simulations indicate that in den- 
sity profiles approximating those of galactic molecular cloud cores 
or cosmological minihalos at high redshift, the I-front remains D- 
type for the lifetime of typical UV sources (Whalen & Norman 
2008a). If rs,o is instead equal to or greater than ro, the I-front 
may briefly convert to D-type, but then rapidly reverts to R-type 
and flash-ionizes the cloud on timescales shorter than the dynami- 
cal time of the gas. Now completely ionized and nearly isothermal, 
strong pressure gradients form wherever there are steep density gra- 
dients, the sharpest of which are found at the edge of what was once 
the edge of the core. These pressure gradients drive the gas outward 
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Figure 6. Test 5 (H II region expansion in an initially-uniform gas): Images of the pressure, cut through the simulation volume at coordinate z = at time 
t = 500 Myr for (left to right and top to bottom) Capreole+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



into the ionized cloud forming a shock that moves with a roughly 
constant velocity in density profiles (Franco et al. 1990). 

In Test 6 we examine the former case, in which the initial 
Stromgren radius is smaller than the core radius. The aim of this 
test is to study the initial transition of the I-front from R-type 
to D-type and back to R-type over a fairly restricted range of 
radii, rather than its long-term behavior thereafter. Accordingly, 
we adopt the following numerical parameters: computational box 
length L = 0.8 kpc, no = 3.2 cm~^, ro = 91.5 pc, zero initial 
ionization fraction, ionizing photon emission rate — 10^° pho- 
tons and initial temperature T — 100 K. The source position is 
at the corner of the computational volume (xs, ys^Zs) — (0, 0, 0). 
Boundary conditions are reflective for the boundaries which con- 
tain the origin and transmissive for the other boundaries. For these 
parameters the I-front changes from R-type to D-type inside the 
core. Once the front reaches the core edge it will accelerate as 
it propagates down the steep density slope. The initial recom- 
bination time inside the core (assuming ionized gas temperature 
T — 10^ K) is tree, core = 0.04 Myr. The ionizing spectrum is 
again that of a 10^ K black body, as expected for a massive, metal- 
free Pop III star. Hydrogen line cooling, recombinational cooling, 
and bremsstrahlung cooling are all included, but again not Comp- 
ton cooling. For simplicity, gravitational forces are ignored and no 



hydrostatic equilibrium is imposed on the cloud. Unlike in Test 5, 
left on their own the pressure forces will accelerate gas outward 
in this density- stratified cloud, albeit those forces are much infe- 
rior than the stronger ones due to pressure from the photoheated 
gas. The running time is tsim = 75 Myr. The required outputs are 
neutral fraction of hydrogen, gas number density, temperature and 
Mach number on the grid at times t = 1,3, 10, 25 and 75 Myr, 
and the I-front position (as defined in Test 5) and velocity vs. time 
along the x-axis. 



3.3 Test 7: Photoevaporation of a dense clump 

In Test 7, a plane-parallel I-front encounters a uniform spheri- 
cal clump in a constant background density field. This problem 
has been studied in many contexts, e.g. in relation to the photoe- 
vaporation of dense clumps in planetary nebulae (Mellema et al. 
1998). Depending on the assumed parameters the clump may ei- 
ther initially trap the I-front, or be flash-ionized without ever trap- 
ping the I-front, the so-called 'cloud-zapping' regime (c.f. Bertoldi 
1989). The condition for an I-front to be trapped by a dense clump 
with number density uh can be derived by defining a "Stromgren 
length", ^s(t), at a given impact parameter r using equations (6) 
and (7), and solving them for each impact parameter (Shapiro et al. 
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Figure 7. Test 5 (H II region expansion in an initially-uniform gas): Images of the temperature, cut through the simulation volume at coordinate z ■ 
t = 500 Myr for (left to right and top to bottom) Capreole+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



: at time 



2004). We can then define the "Stromgren number" for the clump as 
Ls = 2rciump/^s(0), where rdump is the clump radius and is{0) 
is the Stromgren length for zero impact parameter. If > 1, 
then the clump is able to trap the I-front, while if < 1, the I- 
front quickly ionizes the clump and is never trapped. For a uniform 
clump equation (8) reduces to 



is 



(13) 



D-type, in contrast to the less interesting for us "cloud-zapping" 
regime in which the front flash-ionizes the cloud and remains R- 
type throughout. The computational box size is Xbox = 6.6 kpc, 
the radius of the clump is rdump = 0.8 kpc, and its center is at 
(xc, Zc) = (5, 3.3, 3.3) kpc = (97, 64, 64) cells. Hydrogen line 
cooling, recombinational cooling, and bremsstrahlung cooling are 
included, but not Compton cooling. Boundary conditions are trans- 
missive for all grid boundaries. 



and Ls becomes 



2rclumpQf 



(2) 2 



(14) 



The numerical parameters for Test 7 are the same as for Test 
3 in Paper I: a constant ionizing photon flux of F = 10^ s~^cm~^ 
is incident at ^ = 0, the ambient hydrogen gas number density 
and temperature are riout = 2 x 10~^ cm~^ and Tout = 8, 000 
K, respectively, while the initial clump density and temperature are 



^dun 



200 no 



0.04 cm"'' andTduE 



40 K. These pa- 



rameters ensure that outward pressures in the clump balance those 
from the hot gas so that the clump is initially in pressure equi- 
librium with the surrounding medium. The column density of the 
clump is sufficient to trap the I-front and compel its transition to 



With hydrodynamics the evolution beyond the trapping phase 
proceeds very differently from the static Test 3 in Paper I. As 
the heated and ionized gas is evaporated and expands towards the 
source, its recombination rate falls and it attenuates the ionizing 
flux less. As a consequence, the I-front slowly consumes the clump 
until it photoevaporates completely. The required outputs are H I 
fraction, gas pressure, temperature, and Mach number at times 
t — 1, 5, 10, 25 and 50 Myr and the position and velocity of the 
I-front along the axis of symmetry. 
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Figure 8. Test 5 (H II region expansion in an initially-uniform gas): Images of the H II fraction, cut through the simulation volume at coordinate z = at time 
t = 500 Myr for (left to right and top to bottom) Capreole+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



4 RESULTS 
4.1 Test 5 

We start by comparing the fluid flow and ionization structure at two 
characteristic stages of the evolution, at t = 100 Myr, of order of 
one recombination time, which is the start of the I-front conver- 
sion from R- to D-type, shortly before the initial Stromgren radius 
is reached, and at t = 500 Myr, corresponding to a few recombi- 
nation times, when the I-front is D-type preceded by a shock. In 
Figures 2 - 4 and 5 - 7 we show image cuts at coordinate z = 
of the neutral hydrogen fraction, pressure, and temperature at 100 
Myr and at 500 Myr, respectively, while in Figures 8 and 9 we 
show the ionized fraction and number density at 500 Myr. We note 
here that unlike the other simulations which are fully 3-D in both 
the hydrodynamic and the radiative transfer treatment, the RHID 
results are 1-D spherically- symmetric Lagrangian profiles mapped 
onto the 3-D Cartesian grid required in this study. 

With a few exceptions, discussed below, all results exhibit 
reasonably good agreement throughout the flow evolution. As we 
found also for the static tests in Paper I, the majority of the dif- 
ferences are a consequence of the different handling of the energy 
equation and the hard photons with long mean free paths. These 
variations yield different spatial structures in the temperatures (Fig- 



ures 4, 7 and 13) and ionized fractions in the gas just ahead of the 
I-front (which are dictated by the hard photons and non-equilibrium 
chemistry. Figures 8 and 11), but very similar ionization profiles in- 
side the H II region (which sees the whole spectrum of photons and 
is mostly chemically equilibrated. Figures 2, 5 and 11). 

Regardless of the variations in the temperature and ionization 
profiles among the codes, the overall differences in I-front position 
and velocity are very modest, of order only a few percent, with the 
exception of Enzo-RT and, to a lesser extent, HART. The hydro- 
dynamical profiles also cluster fairly closely together. The codes 
basically agree on the temperature structure of the evolving H II 
region over time except for HART, which predicts flat, lower tem- 
peratures at later times, and C^-Ray, which yields higher ionized 
gas temperatures close to the ionizing source due to its simplified 
method for handling the energy, and again Enzo-RT because of its 
monochromatic spectrum. The reason for the sharp drop in pressure 
at 0.6 Lbox at 10 Myr in the HART results is unclear. 

Apart from the differences discussed above, there are several 
features of the HART, LICORICE and Enzo-RT methods worth 
noting. The OTVET moment radiative transfer method used in 
HART is somewhat diffusive, as was already noted in Paper I, 
which results in thicker I-front and less sharp flow features over- 
all. There are some radial striations visible in the LICORICE re- 
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Figure 9. Test 5 (H II region expansion in an initially-uniform gas): Images of the gas number density, cut through the simulation volume at coordinate z - 
at time t = 500 Myr for (left to right and top to bottom) Capreole+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, Flash-HC and Enzo-RT. 



suits, especially in the temperature images that are reminiscent of 
those observed in the CRASH code results in Test 2 of Paper L 
Since LICORICE adopts the Monte Carlo radiative transfer found 
in the original version of CRASH, the radial artifacts in its tem- 
peratures are similarly due to the noise in that version's energy 
sampling scheme, which has been corrected in the latest release 
of the CRASH code (Maselli et al. 2009). The wall effects in the 
upper left and lower right corners of the box in the HART pres- 
sure and Mach number images reflect the fact that mirror rather 
that transmissive boundary conditions were utilized. This is due to 
the natively-periodic nature of the OTVET method, which demands 
special handling in order to run the non-periodic test problems in 
this comparison. The LICORICE, and to a lesser extent the RSPH 
Mach numbers exhibit a somewhat grainy structure deep inside the 
H II region not visible in the other quantities. The origin of these 
features is likely due to the low SPH resolution in the evacuated 
interior of the H II region, which is nearly an order of magnitude 
lower in density than its surroundings. The difference in the degree 
of graininess between the two SPH codes may in part be due to 
how each code's particle data was mapped onto the Cartesian grid. 
The origin of the third outermost band in the RSPH Mach numbers, 
which is not present in those of the other codes, is likely due to the 



utilization of a larger box in that calculation compared to the other 
cases, which changes the flow boundary conditions. 

Several important questions arise in this Test. First, does the 
broadening of the front by high-energy photons from a hard UV 
source alter its radius as a function of time in comparison to a 
monochromatic front with the same average ionized gas temper- 
atures? This issue is key because it determines if the extensive ap- 
proximate analytical solutions to hydrodynamical I-front transport 
that exist in the literature apply to ionization fronts in which there 
is spectral hardening due to hard UV sources. Second, how does 
the penetration of hard photons into the dense shocked neutral gas 
ahead of the I-front alter its structure and flow? Third, how do these 
changes to the shocked flow alter its own rate of advance and that 
of the front? Finally, what are the origin of the distinctive double 
peaks in density, velocity and Mach number in the full spectrum 
I-fronts at intermediate times, and why are they absent in the Enzo 
profiles? 

We cannot resort to comparison of the present code results 
alone to resolve these questions because they are all are multifre- 
quency in nature except for Enzo, and even Enzo integrates over 
the blackbody spectrum to implement the grey approximation to 
radiation transport. These issues can only be settled by comparing 
the multifrequency I-front in Test 5 to a monochromatic one whose 



Cosmological Radiative Transfer Comparison II 




Figure 11. Test 5 (H II region expansion in an initially-uniform gas): Spherically-averaged profiles for ionized fractions x and neutral fractions xhi = 1 — 
at times t = 10 Myr, 200 Myr and 500 Myr vs. dimensionless radius (in units of the box size). 
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Figure 12. Test 5 (H II region expansion in an initially-uniform gas): Spherically-averaged profiles for pressure, p, at times t = \Q Myr, 200 Myr and 500 
Myr vs. dimensionless radius (in units of the box size). 




Figure 13. Test 5 (H II region expansion in an initially-uniform gas): Spherically-averaged profiles for temperature at times t = 10 Myr, 200 Myr and 500 
Myr vs. dimensionless radius (in units of the box size). 




J^/kox J^/kox J^/kox 

Figure 14. Test 5 (H II region expansion in an initially-uniform gas): Spherically-averaged profiles for the hydrogen number density, n, at times t = 10 Myr, 
200 Myr and 500 Myr vs. dimensionless radius (in units of the box size). 
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Figure 16. ZEUS-MP Test 5 ionized fraction (left) and temperature (right) profiles with monoenergetic photons (dashed) and a 10^ K blackbody spectrum 
(solid) at times t = 10 Myr (left pairs), 200 Myr (central pairs) and 500 Myr (right pairs) vs. dimensionless radius (in units of the box size). 



photon energy has been adjusted to yield the same average ionized 
gas temperature as for the 10^ K blackbody spectrum. This guaran- 
tees that any discrepancy in position between the two fronts will be 
due only to the broadening of the front and its modifications to the 
shocked flow just beyond it, not to differences in the average sound 
speed within the H II regions, which is primarily what determines 
the rate of advance of the I-front when it is D-type. This approach 
also ensures that any variations in the structure of the shocked flow 
between the two I-fronts are due to spectral hardening only, since 
both are being driven by the same ionized gas pressure. 

To investigate these points and determine the origin of some 
of the features in the hydrodynamic profiles in Figures 11 -15, we 
performed two fiducial runs of Test 5 with ZEUS-MP. The first 
was with the original 10^ K blackbody spectrum and the second 
was with monoenergetic photons at 17.0 eV. Both had the same 
ionizing photon rate = 5 x 10^^ s~^. The 17.0 eV monochro- 



matic photons establish the same average ionized gas temperature 
as in the multifrequency H II region in ZEUS-MP. We show ion- 
ized fractions, temperatures, velocities, and densities for the two 
runs at t =10, 200, and 500 Myr in Figures 16 and 17. The broad- 
ening of the I-front in the multifrequency calculation is apparent 
at all three times in the ionized fractions, becoming greater as the 
front expands. In contrast, the monoenergetic I-front remains sharp, 
intersecting the multifrequency front at very nearly the same ion- 
ized fraction at all three radii. Except for small differences in the 
elevated values near the source, the two H II regions exhibit nearly 
identical temperatures out to where the full- spectrum I-front broad- 
ens. Out to this same radius the density and velocity profiles are 
also nearly identical. 

The ionization profiles demonstrate that the position of the I- 
front as a function of time is not signficantly altered by either the 
broadening of the front or the partial ionization and heating of the 
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Figure 17. ZEUS-MP Test 5 density (left) and gas velocity (right) profiles with monoenergetic photons (dashed) and a 10^ K blackbody spectrum (solid) at 
times t = Myr (left pairs), 200 Myr (central pairs) and 500 Myr (right pairs) vs. dimensionless radius (in units of the box size). 
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Figure 18. Test 5 (H II region gasdynamic expansion in an initially-uniform 
gas): The evolution of the position and velocity of the I-front. Solid lines 
show the approximate analytical solution as described in the text. Dotted 
horizontal line indicates the approximate value of vr. 



dense shocked gas in front of it by high-energy photons, at least for 
the radii considered in this problem. This affirms that the global dy- 
namics of the D-type I-front depend primarily on the temperature 
(and hence sound speed) of the ionized gas. Past tests by one of us 
of D-type I-fronts in density gradients confirm that this holds 
well beyond the R to D transition surveyed in this Test. Likewise, 



these two models demonstrate that the ionized flow within the H 
II region is also mostly unchanged by spectral hardening over the 
radii enclosed by the computational box. However, here it is impor- 
tant to distinguish between the motion of the I-front and that of the 
shocked flow it drives. The latter is dramatically altered by spectral 
hardening as we discuss below. 

The simple- structured 3000 K layer of shocked gas driven by 
the monochromatic I-front is split into the double-peaked structure 
in the full- spectrum front at t = 200 Myr and is present in the pro- 
files of all the codes except Enzo, as shown in Figures 14 and 15. 
This feature is transient and disappears by t = 500 Myr. Its origin 
is the heating of the dense shell by the high-frequency photons. At 
200 Myr, they partially ionize the base of the shocked shell: ion- 
ized fractions of 10% or more extend out to 0.5 Lhox- The high fre- 
quency tail of the spectrum cannot maintain large ionized fractions 
in this layer but does effectively deposit heat there, as evidenced 
by the rise in temperature at 0.45 Lhox, which is positioned ap- 
proximately in the valley between the two peaks in the density and 
velocity profiles. This energy ablates the lower layer of the dense 
shocked shell, driving both inward and outward photoevaporative 
flows in the frame of the shock that split the density and velocity 
peaks into two smaller ones. The forward flow accelerates the gas 
in the outer peak to 7 km s~^ by 200 Myr. However, pressure gradi- 
ents from the ionized interior of the H II region drive the inner peak 
to higher velocities that cause it to later overtake the forward peak 
(500 Myr). At this distance from the central source, high energy 
photons do heat the base of the shocked shell but not to sufficient 
temperatures to create backflow, as seen in the disappearance of the 
temperature bump that was present at 200 Myr. However, they do 
smear out the sharp interface between the ionized and shocked gas 
temperatures that is present in the monoenergetic front at 500 Myr. 
In contrast, monoenergetic photons result in much simpler structure 
in both the front and the dense shell at 200 and 500 Myr. At 500 
Myr, ~ 15,000 K ionized gas drives a clearly defined shocked shell 
and there are no ablation flows. The absence of backflows allow 
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Figure 19. Test 6 (H II region gasdynamic expansion down a power-law 
initial density profile): The evolution of the position and velocity of the I- 
front. 

peak gas velocities to reach higher values in the shell at intermedi- 
ate times than in the hard spectrum case. 

What effect does pre-heating by hard photons leaking ahead of 
the I-front have on the propagation of the shocked flow? It weakens 
the shock, as evidenced by the smaller density jump, lowering the 
density compression there and thus enlarging its detachment from 
the I-front. This can be clearly seen by comparing the position of 
the shock for the two I-fronts in Figure 17. Thus, while the positions 
of the two fronts are nearly identical, the shocked flow of the mul- 
tifrequency front is well ahead of that of the monochromatic one. It 
is clear from this comparison that multifrequency photon transport, 
or the use of lookup tables of ionizing rates as a function of opti- 
cal depth as a proxy, is necessary to capture the correct structure of 
I-fronts and shocks driven by high-temperature UV sources. 

Having isolated the effect of spectral hardening on both the 
dynamics of the I-front and on the shocked neutral flows it drives, 
we can now disentangle the sometimes competing effects that ac- 
count for the differences that do exist among the hydrodynamic 
profiles in Test 5, both among the codes with multifrequency 
physics and between those codes and Enzo. Setting aside the Enzo 
results for the moment, the minor spread in I-front position in the 
other codes can now be traced to the variations of their H II region 
temperatures. This is primarily due to how each code handles the 
energy equation. Furthermore, it is now clear that (1) the origin of 
the temperature bumps in the shocked gas at intermediate times; 
(2) the double peaks in the densities, velocities, and Mach num- 
bers; and (3) shocks that are more fully detached from the I-front 
are all a consequence of spectral hardening. 

The profiles that are most distinct from the rest are those of 
Enzo-RT. This can now be understood to be due to the lack of 
spectral hardening in this code. Although it carries out an integra- 
tion over the 10^ K spectrum to compute the photoionization cross 
sections, Enzo-RT it employs a grey approximation for the pho- 
toionization cross-section, i.e. a cross-section independent of the 



frequency and is therefore essentially monochromatic in its current 
form. No hard photons means no pre-heating ahead of the front and 
also a much sharper I-front. This lack of pre-heating in turn means 
that, as in the ZEUS -MP monochromatic results above, the Enzo- 
RT shock is stronger than the others, as evidenced by its higher 
Mach number, density compression and pressure jump. Because the 
Enzo I-front remains sharp, its shocked flow also exhibits the single 
peak associated with monoenergetic I-fronts. Because it is stronger, 
the shock in the Enzo profiles propagates somewhat more slowly, 
lagging behind those of the other codes. On the other hand, the 
Enzo-RT I-front actually leads the others, almost coinciding with 
the shock. This is a consequence of the greater average ionized gas 
temperatures in its H II region in comparison to the others, a re- 
sult of its integration over the blackbody spectrum. The Enzo-RT 
temperature profile is similar to the flat profile found by HART, 
but at a higher value. Its origin is unknown but could be associated 
with its flux-limited diffusion radiative transfer, which shares some 
similarities with the OTVET method in HART. 

Except for the variations noted above, the codes agree well on 
I-front position and velocity but are at variance with the analytical 
solution, leading it by roughly 30%. This is in part because the so- 
lution plotted in Figure 17 is for fixed temperatures of 15,000 K 
behind and 1000 K ahead of the I-front (allowing for some pre- 
heating by hard photons, which alters the values of vr and vd), 
while in reality the profiles exhibit a complex temperature struc- 
ture. Furthermore, as we discussed above, the analytical solution 
describes well the early and late evolution, but not the intermediate 
one. The H II region radial evolution at late times does not exactly 
match the asymptotic t^^^ slope predicted for self-similar flows but 
approaches it at larger radii. This is to be expected since the box 
size was chosen to enclose only the transition of the I-front from 
R-type to D-type, during which the assumption of self- similarity is 
not satisfied. The codes asymptotically approach the expected so- 
lution as the fronts grow in radius. 

Finally, we note that the I-fronts in this test are dynamically 
stable. If H2 cooling, LW photodissociation, and self- shielding to 
LW photons had been included, violent hydrodynamical instabili- 
ties mediated by H2 cooling might have erupted in the fronts after 
becoming D-type, as explained in greater detail in Test 6 below. 
Line cooling in H alone appears to be unable to incite such insta- 
bilities (Whalen & Norman 2008b). 

4.2 Test 6 

We start our analysis with a head-to-head comparison of the evolu- 
tion of the position and velocity of the I-front, plotted in Figure 19. 
In the velocity plot we clearly see the evolution stages outlined in 
the discussion of Test 6 above. Initially, while it is still within the 
density core the I-front moves very fast (is of R-type), but precip- 
itously slows down as it approaches its Stromgren radius (whose 
precise value is temperature-dependent, but is slightly smaller than 
the core radius chosen here). The fast R-type phase is over within 
a fraction of a Myr, after which the expansion becomes pressure- 
driven, and the front itself converts to a D-type led by a shock. The 
I-front speed reaches a minimum of just a few kilometers per sec- 
ond, well below vr - the critical velocity defined in § 3.1. We note 
that although some of the results appear to never show I-front ve- 
locities below ^10 km/s, this is in fact due to insufficient number 
of early-time snapshots being saved in the I-front evolution data. 
For this reason the short transition stage does not appear in some 
of the plotted results, and this does not imply any problem with the 
codes. The later- time evolution is not affected, as long as the ac- 
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Figure 20. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Images of the H I fraction, cut through the simulation volume 
at coordinate z = at time t = 25 Myr for (left to right and top to bottom) Capreole+C^ -Ray, TVD-^C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, 
and Flash-HC. 



tual time- stepping in the code is sufficiently fine to properly follow 
the early evolution. Once out of the core, the I-front re- accelerates 
as it descends the steep density gradient, eventually reaching 
speeds of 25 — 28 km/s. These speeds never surpass vr and there- 
fore the front remains D-type until leaving the computational vol- 
ume at t 30 Myr 

All codes agree on the later-time, pressure-driven expansion, 
both qualitatively (a slow, D-type I-front, preceded by a relatively 
weak shock, as we shall see below), and quantitatively. There are 
some modest differences in the I-front speed, ^ 10% or less be- 
tween cases, which results in I-front positions whose spread grows 
with time, but never exceeds ^ 5 — 7%. 

Next we turn our attention to the overall structure of the fluid 
flow and ionization, shown in 2-D cuts along the x-y plane of the 
H I and H II fractions, density, temperature and Mach number at 
time t — 25 Myr in Figures 20-24. We again remind the reader that 
unlike the other simulations which are fully 3-D in both the hydro- 
dynamic and the radiative transfer treatment, the RHID results are 
1-D spherically- symmetric Lagrangian profiles mapped onto the 3- 
D Cartesian grid required in this study. There is good agreement 
between the results, in terms of the positions of the I-front and the 
shock, the size of the growing H II region and its ionization, density 



and temperature structures. There are some differences in the level 
of hard photon penetration ahead of the I-front and the temperature 
distribution, similar to the ones we observed in Test 5 and Paper I. 

However, in this test also a new kind of difference shows up, 
namely the appearance of instabilities near the I-front. Such insta- 
bilities occur for several of the codes, and their nature varies be- 
tween codes. In the cases of C^-Ray-i-Capreole and LICORICE, 
the instabilities are clearly visible in the ionized fractions, tem- 
peratures, densities, and Mach numbers, while in Flash-HC they 
are mostly visible in the temperatures and ionized fractions. RSPH 
exhibits a minor anomaly in only the temperature at 25 Myr. The 
RHID data cannot exhibit such instabilities because they are ID 
spherical polar coordinate profiles mapped onto the 3D Cartesian 
grid mandated for this test. The ZEUS-MP profiles, which were 
computed on a 3D spherical polar coordinate grid and then mapped 
onto 3D Cartesian coordinates, manifest no instabilities in any of 
the profiles, and the HART results do not show them either. Are 
these instabilities physical or numerical? 

Three types of dynamical instabilities in ionization fronts have 
been discovered in the past thirty years. The first type occurs in D- 
type ionization fronts whose shocked neutral gas shells can cool 
efficiently by radiation (Giuliani 1979; Garcia- Segura & Franco 
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Figure 21. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Images of the H II fraction, cut through the simulation volume 
at coordinate z = at time t = 25 Myr for (left to right) and top to bottom) Capreole+C^-Ray, TVD+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, 
and Flash-HC. 



1996). Cooling collapses the gas into a cold thin dense layer that is 
prone to oscillations and fragmentation (Vishniac 1983). Ionizing 
UV radiation then opportunistically escapes through the cracks in 
the shell and flares outward in violent instabilities. However, in the 
current Test only H lines can cool the shell, and recent numerical 
experiments prove that such cooling is too inefficient to initiate dy- 
namical instabilities of this type (Whalen & Norman 2008a). The 
fact that thin- shell instabilities do not arise in the Test 5 profiles 
further attests to the fact that H line cooling is not responsible for 
the corrugations in the Test 6 profiles. In general, shocks that accel- 
erate down power law density gradients steeper than r~^ are also 
prone to Rayleigh-Taylor instabilities that do not require radiative 
cooling. They too are also capable of inciting violent instabilities in 
I-fronts, but are not relevant to the r~^ gradients in the current test. 
Another type of instability can appear in D-type fronts when pho- 
tons are incident to the front at oblique angles (Williams 2002), but 
they cannot develop in I-fronts given the imposed initial spherical 
symmetry of this test. 

The third type, shadow instabilities, can appear when a density 
perturbation is advected through an R-type front, forming dimples 
that erupt into violent instabilities when the I-front becomes D-type 
(Williams 1999). Although the density profile in this test is radially 



symmetric, prescriptions for imposing spherically- symmetric pro- 
files on a Cartesian mesh as a function of radius that are too simple 
can lead to minor departures from radial symmetry in densities be- 
tween neighbor grid points. As the front crosses these mesh points 
it can become dimpled just as if real physical perturbations had tra- 
versed it. These corrugations would then grow into the much more 
prominent features visible in the 25 Myr images upon transforma- 
tion of the front to D-type. If this were the case, they would be 
dampened by employing higher grid resolution or by a better pre- 
scription for smoothing densities between neighbor mesh points. 
We tested the latter possibility by applying an algorithm in C^-Ray 
that sub-sampled volumes enclosed between adjacent points in ra- 
dius with thirty more finely subdivided shells and then interpolated 
densities accordingly between the points. This measure failed to al- 
leviate the instabilities in the I-front in C^-Ray, suggesting that they 
are not shadow instabilities. 

Instead, their shape and position at early times reveal that 
they are the infamous 'carbuncle' phenomenon or 'odd-even de- 
coupling' (Quirk 1994). Low-diffusion solvers, such as the Roe's 
approximate Riemann solver employed in Capreole-i-C^-Ray and 
the PPM Riemann solvers in Flash-HC and Enzo-RT, sometimes 
suffer from these numerical instabilities. They occur when shocks 
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Figure 22. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Images of the density, cut through the simulation volume at 
coordinate z = at time t = 25 Myr for (left to right and top to bottom) Capreole+C^ -Ray, TVD+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, 
and Flash-HC. 



travel parallel to one of the coordinate axes, in this case beginning 
near the symmetry axes of the expanding shell. These instabilities 
are well-known and mostly understood, although their occurrence 
is not always predictable. Although by 25 Myr most of the shell 
has been disrupted, we find that the perturbations begin near the 
axis and exhibit the characteristic morphology of the carbuncle in- 
stability. While the usual solution is to artificially introduce more 
diffusion only where it is needed, this approach did not suppress 
the phenomenon in C^-Ray. If the test is instead run with C^-Ray 
coupled to the TVD solver of Trac and Pen (Trac & Pen 2004), 
which is more diffusive and not known to suffer from the carbuncle 
instability, the shell indeed remains well behaved. We also note that 
this carbuncle instability, while fairly violent, does not in fact affect 
the fluid flow or I-front propagation significantly and the I-front po- 
sition and velocity evolution discussed above and the spherically- 
averaged profiles discussed below are still in good agreement with 
the other results. Therefore, effect of this is fairly modest at the 
early times studied here, but as the previously cited work demon- 
strates, they could become important at later times. 

Another possibility is that there is a hitherto unknown break- 
out instability associated with the transition of the D-type I-front 
back to R-type as it descends the density gradient. However, such 



an instability would not have the opportunity to propagate through- 
out and disrupt the entire shell during breakout as in the C^-Ray 
results because the transition from D-type to R-type is too abrupt, 
and, as discussed above, the I-front remains D-type within the com- 
putational volume. Such an instability would instead be manifest 
as a premature supersonic runaway of radiation from the surface 
of the shell along certain lines of sight with little disruption of the 
shell itself, as observed in the Flash-HC results. However, this does 
not happen in the ZEUS-MP profiles, which are computed on a 
3D spherical coordinate grid that is naturally suited to spherically- 
symmetric density fields and on which the 'corner' effects inher- 
ent in Cartesian grids are absent. Furthermore, as the spherically 
averaged profiles show, the I-front in this test is not on the verge 
of breaking past the shock and becoming R-type at late times. We 
therefore conclude that the instability in the Flash-HC results is not 
physical and that the early breakout of radiation there is at least 
partly the result of gridding a spherical density on a Cartesian grid. 
The effect of this is fairly modest, however, and does not disturb 
the overall dynamics significantly. 

The irregular morphology of the shell in the LICORICE re- 
sults and to a small extent in the RSPH profile at 25 Myr is likely 
due to spurious fluctuations in the density field where the local 
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Figure 23. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Images of the Mach number, cut through the simulation volume 
at coordinate z = at time t = 25 Myr for (left to right) and top to bottom) Capreole+C^ -Ray, TVD+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, 
and Flash-HC. 



particle number changes sharply, in this case in the vicinity of 
the dense shell. This well-known feature of SPH, as discussed in 
section 2.7, is what probably allows radiation to preferentially ad- 
vance along lines of sight through low-density fluctuations in the 
two profiles. The larger effects for LICORICE compared to RSPH 
are probably due to the usage of a grid to perform the radiative 
transfer in the former. Once again, none of these effects appears to 
affect the overall evolution significantly, but they might matter in 
certain astrophysical situations. 

The HART results exhibit banding in all the profiles except for 
HI fraction at 25 Myr. The origin of these features is unclear, but is 
possibly related to the much coarser AMR griding used around the 
outer edges. They may also be related to the greater diffusivity of 
the OTVET algorithm, although no such features were observed in 
the other tests performed with OTVET. However, the density profile 
found by HART is much flatter, with no clear dense shell swept by 
the shock, in clear contrast to all other results. It is possible that the 
time step applied to the gas energy updates in HART is too coarse 
for I-fronts in r~^ density gradients, which has been found to lead 
to banding in temperatures and densities in H II regions in stratified 
media (Tenorio-Tagle et al. 1986; Whalen & Norman 2006). 

We finally note that had H2 cooling been included in either 



this test or in Test 5, violent physical instabilities might have arisen 
in the I-front when it became D-type (Whalen & Norman 2008a). 
Hard UV spectra significantly broaden ionization fronts, forming 
regions of a few thousand K and ionized fractions of 10% in their 
outer layers. These are prime conditions for the catalysis of H2 via 
the H~ channel, which forms between the I-front and the dense 
shocked shell when the front becomes D-type (Ricotti et al. 2002). 
H2 - H, H2 - e~ , and H2 - H^ collision channels emit ro-vibrational 
lines (Lepp & Shull 1983; Galli & Palla 1998b; Glover & Abel 
2008) that can radiatively cool the base of the shocked layer and 
incite dynamical instabilities there just as metal ions do in galactic 
I-fronts. These instabilities may have been common in the early 
universe, such as in UV breakout from the first star-forming clouds. 
They appear if there is enough H2 in the cloud to self- shield from 
the Lyman- Werner (LW) flux (11.18-13.6 eV) also being emitted 
by the source, which photodissociates molecular hydrogen (Draine 
& Bertoldi 1996). Thus, while the instabilities manifested by some 
of the codes in Test 6 are numerical, physical ones are possible 
when H2 cooling, LW photons, and self- shielding to LW radiation 
are properly included. Since not all the codes contain these physical 
processes, it was not included in any of the current tests, but will be 
a target for future stages of this comparison project. 
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Figure 24. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Images of the temperature, cut through the simulation volume 
at coordinate z = at time t = 25 Myr for (left to right and top to bottom) Capreole+C^ -Ray, TVD+C^-Ray, HART, RSPH, ZEUS-MP, RHID, LICORICE, 
and Flash-HC. 
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Figure 25. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Spherically-averaged profiles for ionized fractions x and 
neutral fractions xhi = I — x at times t = 3 Myr, 10 Myr and 25 Myr vs. dimensionless radius (in units of the box size). 
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Figure 26. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Spherically-averaged profiles for the gas number density, n, at 
times t = 3 Myr, 10 Myr and 25 Myr vs. dimensionless radius (in units of the box size). 




Figure 28. Test 6 (H II region gasdynamic expansion down a power-law initial density profile): Spherically-averaged profiles for temperature at times t = 
3 Myr, 10 Myr and 25 Myr vs. dimensionless radius (in units of the box size). 
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Figure 30. Test 7 (Photoevaporation of a dense clump): The evolution of 
the position and velocity of the I-front along the axis of symmetry through 
the centre of the clump. 



In spite of the prominence of the numerical instabilities in 
some of the codes, we re-iterate that they were not catastrophic 
to the overall dynamics over the range of radii and interval of time 
which we study here, as shown by the spherically averaged hydro- 
dynamical profiles. In Figs 25 - 29 we show ionization fractions, 
number densities, pressures, temperatures and Mach numbers at 3, 
10, and 25 Myr. Comparison of ionization fractions and densities 
at 3 Myr indicates that the I-front is D-type at r 120 pc, some- 
what beyond the flat central core of the initial density profile. A thin 
layer of shocked neutral gas is visible at r ^ 140 pc in the Mach 
number profile. Acoustic waves are evident at r < 100 pc within 
the H II region in both the density and Mach number plots and are 
consistently reproduced by all the codes. 

At this early stage, the I-front widths (customarily defined by 



the difference between the positions at which 0.1 and 0.9 ionization 
fractions are reached) vary from ^ 20 pc to 40 pc. At the relatively 
high inner-profile density of n ^1 cm~^, the mean free paths (mfp) 
of 13.6 eV and 60 eV photons, which roughly bracket the avail- 
able energies in the 10^ K black-body spectrum used for this test, 
are ^ 0.05 pc and 4 pc, respectively. The intrinsic width of the I- 
front is approximately 20 mean-free paths, or between 1 and 80 pc. 
Therefore, all the codes give widths roughly consistent with the ex- 
pected values, but somewhat on the wider side, primarily due to the 
diffusivity of some of the algorithms. In particular, LICORICE and 
HART have wider fronts, while ZEUS-MP has the narrowest one 
and the rest are spread between those two extremes. As explained 
above, the details of the structure of the I-front are interesting since 
they relate to the formation of molecular hydrogen in its outer lay- 
ers (e.g. Ricotti et al. 2002) where hard, deeply penetrating photons 
could yield a positive feedback mechanism during early structure 
formation. The post-front (ionized) gas temperatures of the codes 
at these early times differ by at most 10%. The low Mach num- 
bers outside the H II region at early times, ranging from 0.001 to 
0.01 arise because hydrostatic equilibrium was not imposed on the 
original density profile in any of the codes except for ZEUS-MP. 
Pressure forces gently accelerate the gas outward, but this has little 
effect on the late-time evolution of the H II region. For the purposes 
of this comparison the lack of initial hydrostatic equilibrium is ir- 
relevant, as long as all codes start from the same initial conditions. 

By 10 Myr the H II region has grown to 240 pc, with all results 
agreeing well on the I-front position. There is very little variation in 
the ionization structure inside the H II region, with only LICORICE 
finding a slightly lower level of ionization. The differences in the 
pre-front ionization structures are much more pronounced, underly- 
ing again the variety in the treatments of multi-frequency photons. 
All the codes still find postfront gas temperatures within 10% of 
one another at most radii. The temperature profiles drop sharply 
just beyond the I-front as before, but then briefly plateau at 10^ K 
for ~ 16 pc before falling further. This is the dense ambient neutral 
gas shell- swept up by the shock, clearly seen in the density profiles 
(Figure 26), which is sufficiently hot and dense to become colli- 
sionally ionized to a small degree. However, the minute residual 
ionized fractions (10~^ - 10~^) in and beyond the shell in most 
of the plots occur because the I-front broadens over time. As more 
neutral gas accumulates on the shell with the expansion of the H II 
region, its density decreases because its area grows, and its opti- 
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cal depth to photons in the high-energy tail of the spectrum de- 
creases. ZEUS-MP still finds the sharpest I-front and LICORICE 
the thickest, with the rest dispersed between. No single cause can 
be ascribed to the moderate variation in I-front structure amongst 
the codes; for example, both ZEUS-MP and RHID perform mul- 
tifrequency ray-tracing radiative transfer with similar integration 
schemes and frequency binning, but RHID has a noticeably wider 
I-front. Tabulating pre-computed frequency-dependent ionization 
rate integrals as a function of optical depth as an alternative to full 
multifrequency RT in Capreole-i-C^-Ray leads to a somewhat dif- 
ferent structure for the front. There is an unmistakable trend toward 
greater diffusivity with the SPH-coupled radiative transfer codes 
RSPH and LICORICE that is likely related to the inherent difficulty 
in representing low-density regions with SPH particles and the ten- 
dency of SPH to broaden shocks. Nonetheless, the grid-based codes 
and RSPH agree to within a few percent on the density structure 
of the shocked shell at 10 Myr. LICORICE does not resolve the 
shell as well, but this would likely be remedied by using more par- 
ticles to resolve the flow, or by using a more adaptive smoothing 
kernel. Overall, the codes agree reasonably well on the shock po- 
sition and the corresponding density and pressure jumps (although, 
as was discussed above, HART yields higher and more uniform 
density and pressure distributions behind the shock than the other 
the methods, which agree on that quite well among themselves). 

At 25 Myr the I-front is at 640 pc, approaching the boundary 
of our computational volume. At this stage the subsonic expansion 
of the front with respect to the sound speed in the ionized gas is 
evident: acoustic waves have erased density fluctuations up to the 
shocked shell in the pressure and density plots. The acceleration of 
the shock down the density gradient can be seen in the heating by 
the shock: the temperature of the dense shell is 25% greater than at 
10 Myr. The velocities beyond the shock are now 20% of the sound 
speed of the neutral gas, and the peak density of the shell has fallen 
to 0.4 cm~^. There is a 10% variation in the position of the I-front 
among the codes that is not attributable to differences in chemistry 
or radiative transfer because of the uniformity in ionized gas tem- 
peratures (and therefore sound speeds). More likely, it is due to the 
variety of hydrodynamics schemes (both grid and particle based) 
applied to the models. Apart from the diffusivity of some of the al- 
gorithms as manifest in I-front structure, we find good agreement 
on the evolution of the H II region in this stratified medium between 
all the codes. Direct multifrequency RT and approximations to full 
multifrequency transfer with precomputed ionization integrals both 
yield extended I-front structures, but it is difficult to assess which 
is more accurate since even the two direct methods disagree with 
each other to some degree. The disagreement between the multifre- 
quency codes on the width of the I-front is probably due to their dis- 
cretization of the blackbody curve and resultant binning of ionizing 
photon rates according to energy, since they otherwise employ the 
same ionization cross sections and photon-conserving ray tracing. 
The number of bins per decade in energy and their distribution in 
frequency can lead to different thicknesses for the front. They have 
a much smaller effect on the temperature of the ionized gas, and 
therefore the position of the front, because the temperature is more 
strongly governed by the cooling rates than by minor discrepancies 
in spectral profile. 

4.3 Test? 

The evolution of the I-front along the axis of symmetry through 
the center of the dense clump is shown in Figure 30. The I-front 
starts off very fast, R-type in the low-density medium surrounding 



the dense clump, but slows down quickly once it enters the high- 
density gas, which occurs in less than a Myr. Thereafter, the front 
slows down more gradually as source photons encounter more and 
more recombining atoms in the photoevaporative flow, which atten- 
uates the flux which reaches the I-front. This initial trapping phase 
is largely over by t = 1 Myr, yielding a thin ionized layer in the 
dense clump on the source side and a clear shadow behind, as illus- 
trated in Figure 31. Due to the short evolution time, by this point 
the gas is still essentially static and Test 6 reproduces the analo- 
gous stage in Test 3 in Paper I. There are only a few modest dif- 
ferences between the neutral gas distributions. The boundary of the 
shadowed region "flares" especially for the RSPH result, primar- 
ily because their particle neighbour list based ray-tracing scheme 
inevitably introduces some "diffusion of optical depth", whereby 
high optical depth spreads through the neighbour list. Flaring of 
the boundary is also in part due to the interpolation procedures to 
set up the initial conditions which, with their sharp boundaries, are 
unnatural for SPH and thus difficult to represent well. In fact, even 
grid-based codes exhibit similar problems, since the low resolution 
required in this Test imposes some grid artefacts on the spherical 
dense clump that are later manifest as ripples in the ablation shock 
of the clump unless a smoothing procedure is applied at the setup 
of the problem. To minimize artificial features in the photoevapora- 
tive flow from the clump, ZEUS-MP and C^-Ray implemented the 
same smoothing procedure to the clump as in Test 6 in their initial 
density profiles. There is also some faint striping of the neutral frac- 
tion in the low-density region for the case of LICORICE, probably 
due to insufficient Monte-Carlo sampling, as was discussed earlier. 
However, this does not have any apparent effect on the photoevap- 
oration of the clump. 

Once the front speed drops below vr a shock starts to form 
ahead of it, converting it to D-type front. The photoheated mate- 
rial on the source side starts photoevaporating, by blowing super- 
sonic wind towards the ionizing source. The I-front slowly eats its 
way into the dense clump, as shell after shell of gas boils off and 
joins the wind. The I-front velocity gradually drops to a few km/s 
and its position remains roughly constant. Some differences among 
the derived I-front evolution in position and velocity are observed, 
but they remain small throughout the evolution, never exceeding 
10% in terms of position. In Figure 32 we show images of the neu- 
tral hydrogen fraction at t = 10 Myr. Overall results agree fairly 
well, with the expanding wind and shadow at very similar stages of 
evolution. There are also a few, relatively minor, differences which 
should be noted. The RSPH result remains somewhat more diffuse 
and asymmetric then the rest, as noted above, but as the evolution 
proceeds the differences are somewhat less notable. However, the 
photoevaporation does proceed somewhat more rapidly in this case 
due to the inevitably more diffuse initial conditions. There is some 
leaking of light at the edges of the LICORICE shadow which is not 
seen in the other results and should therefore be related to the ra- 
diative transfer method employed, rather than to other factors, e.g. 
to the limb column density of the clump being small and allowing 
some light to go through. Finally, there are some uneven features at 
the edge of the shadow on the source side in the case of Flash-HC, 
whose origin is currently unclear. 

The pressure images at t = 10 Myr in Figure 33 essentially 
agree, with only minor morphological differences between the re- 
sults. The shadow is somewhat thicker and less squeezed at the 
edges for C^-Ray, ZEUS-MP and Flash-HC, compared to RSPH 
and LICORICE, with Coral results intermediate between the two 
groups. The reason for this difference becomes apparent from the 
corresponding temperature images (Figure 34). In the cases of C^- 
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Figure 31. Test 7 (Photoevaporation of a dense clump): Images of the H I fraction, cut through the simulation volume at coordinate z = at time t = I Myr 
for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 




Figure 32. Test 7 (Photoevaporation of a dense clump): Images of the H I fraction, cut through the simulation volume at coordinate z = at time t = 10 Myr 
for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 



Ray, ZEUS-MP and Flash-HC there is clear temperature gradi- 
ent from the edges of the shadow going inward, while for RSPH 
and LICORICE this temperature gradient is largely absent. There 
are also noticeable temperature variations within the clump for 
C^-Ray, ZEUS-MP and Flash-HC which are less pronounced for 
RSPH and LICORICE. The reason for these differences appears to 
be the different levels of penetration by hard photons through the 



high column density material in the clump, and the corresponding 
varying levels of energy deposit by those photons. The variations in 
evolution this introduces seem minor in our particular test problem, 
but such differences might matter more in problems in which the 
precise level of the number of free electrons and the local temper- 
ature within dense clumps is of importance. One example may be 
the study of the the production of molecular hydrogen within dense 
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Figure 33. Test 7 (Photoevaporation of a dense clump): Images of the gas pressure, cut through the simulation volume at coordinate z = at time t = 10 
Myr for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 
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Figure 34. Test 7 (Photoevaporation of a dense clump): Images of the gas temperature, cut through the simulation volume at coordinate z = at time t = 10 
Myr for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 



regions irradiated by UV radiation, which can regulate (stimulate 
or suppress) local star formation (Iliev et al. 2006a; Whalen et al. 
2008a). Finally, the Mach number images at t = 10 Myr are shown 
in Figure 35. The supersonic wind which starts to blow towards 
the ionizing source is clearly visible, with only small differences 
in terms of the thickness of this layer and the Mach number values 
between the different runs. The only peculiarity visible here is that 



in the case of ZEUS-MP this supersonic layer is almost spherical, 
surrounding the clump from all sides, which is not seen in any of 
the other results. This appears to be a consequence of the very cold 
(T < 1000 K) region remaining at the back of the clump, which is 
not present in any of the other cases (see also Figure 43). The rea- 
son for this region remaining so cold in the ZEUS-MP simulation 
is unclear at present, considering that (as we discussed above) the 
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Figure 35. Test 7 (Photoevaporation of a dense clump): Images of the flow Mach number, cut through the simulation volume at coordinate z = at time 
t = 10 Myr for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 



spectrum hardening and penetration of hard photons through the 
clump and into the shadow are similar to C^-Ray, Flash-HC and 
Coral and stronger than RSPH and LICORICE. 

By t = 50 Myr (Figures 36-39) the photoevaporation process 
is well advanced. The region swept by the expanding supersonic 
wind has grown quite large and takes up a significant fraction of 
the simulation volume. There are only modest differences in its size 
between the different codes. In the case of Flash-HC and, to a lesser 
extent Coral, the edge of the expanding wind region is uneven, as 
a consequence of the grid effects in the initial conditions, as dis- 
cussed above, when representing a spherical object on a relatively 
coarse rectangular grid with no interpolation used. These grid ef- 
fects could be seen at earlier times as well, but at a lower level. 
The overall size of the wind region is the same, however, thus this 
problem does not affect the evolution significantly. 

A small core region from the initial clump remains neutral and 
still casts a clear shadow which also remains neutral in all cases. 
This neutral region is moderately compressed by the higher exter- 
nal pressure of the ionized and heated gas surrounding it. The size 
of this neutral dense core and its shadow varies between the runs, 
being somewhat larger for Capreole-i-C^-Ray, Flash-HC and Coral 
than for RSPH, ZEUS-MP and LICORICE. There is also some 
'flaring' (i.e. widening) of the shadow for Flash-HC and Coral, 
probably due to the specific interpolation weighting used in the 
short characteristics methods they both employ (for discussion and 
testing of this see Mellema et al. 2006b, Appendix A). 

In the pressure and temperature images shown in Figures 37 
and 38 we clearly see the shocked shell of gas swept up by the 
supersonic wind of evaporating clump material. The inner zone on 
the side of the clump facing the source is being evacuated and is ac- 
cordingly colder due to adiabatic cooling, while the outer shocked 
shell is much hotter, with temperatures reaching 40,000-70,000 K 
(note the different upper limits for the temperature images). Some 
quantitative and morphological differences are easily noticed. The 



evacuated region yields a shell of low pressure whose depth varies 
between the runs by about an order of magnitude, from the very low 
pressure ^ 10~^^g/cm/s^ found by Flash-HC, through the inter- 
mediate cases of RSPH and ZEUS-MP, to the relatively higher pres- 
sure - 10"^^g/cm/s^ found by Capreole-hC^-Ray, LICORICE 
and Coral. The dense, high-pressure central region which remains 
neutral and the shadow behind it show quite different morpholo- 
gies between the runs, clearly seen in the pressure images (less so 
in the temperature ones, due to the lack of color dynamic range). 
This morphology arises as a consequence of successive reflecting 
oblique shocks which form behind the evaporating clump by the 
interaction between the evaporative wind and the partly collapsed 
shadow squeezed inward by the high external pressure of the ion- 
ized region. The reason for the morphological differences between 
the cases is most probably a slight difference in the timing of these 
shocks for each run, but ascertaining this will require more detailed 
analysis of the evolution. 

Finally, the Mach number images at t = 50 Myr shown in 
Figure 39 show that while the wind is clearly supersonic, with 
Mach numbers of a few, the shocked swept material moves sub- 
sonically. The peak Mach numbers vary from 2 to 3.7, with typical 
peak values around 3. The shock is clearly somewhat weaker for 
Coral, a consequence of this code's more diffusive hydrodynamic 
solver (based on van Leer flux splitting). All other methods, both 
Eulerian grid-based (Capreole-hC^ -Ray, ZEUS-MP and Flash-HC) 
or particle-based (RSPH, LICORICE) yield very similar results in 
terms of Mach number values. The only significant difference be- 
tween the results is again the more spherical high Mach number 
shell found by ZEUS-MP. 

Next we turn our attention to the statistical distributions of 
the temperature (shown in Figure 40) and the Mach number (in 
Figure 41). We notice that three distinct temperature phases, rep- 
resented by the three peaks of the histograms, exist throughout the 
evolution - hot, photoionized gas with temperatures T ^ 25, 000 — 
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Figure 36. Test 7 (Photoevaporation of a dense clump): Images of the H I fraction, cut through the simulation volume at coordinate z = at time t = 50 Myr 
for (left to right and top to bottom) Capreole+C^-Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 




Figure 37. Test 7 (Photoevaporation of a dense clump): Images of the pressure, cut through the simulation volume at coordinate 2; = at time t = 50 Myr 
for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 



45, 000 K, very hot, T > 50, 000 K, shock-heated gas and a cold 
phase, consisting in part of self- shielded gas and in part of adiabati- 
cally cooled gas behind the expanding supersonic wind. These three 
phases are observed in all cases and the histograms are very similar. 
The Mach number histograms are in good agreement as well. For 
RSPH and LICORICE the hot, shocked phase is less distinct from 
the photoionized phase. The shocked gas temperature is a bit higher 



for LICORICE, due to the stronger shock (evidenced by the higher 
peak Mach number) observed in this case. On the other hand, the 
temperatures found by Capreole-i-C^-Ray are somewhat lower than 
the rest, which is related to the more approximate treatment of the 
energy equation in that case. 

Finally, in Figures 42-44 we present cuts along the x-axis of 
the neutral fraction, xhi, temperature, T, pressure, p, and Mach 
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Figure 38. Test 7 (Photoevaporation of a dense clump): Images of the temperature, cut through the simulation volume at coordinate z = at time t = 50 Myr 
for (left to right and top to bottom) Capreole+C^-Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 






Figure 39. Test 7 (Photoevaporation of a dense clump): Images of the flow Mach number, cut through the simulation volume at coordinate z = at time 
t = 50 Myr for (left to right and top to bottom) Capreole+C^ -Ray, RSPH, ZEUS-MP, LICORICE, Flash-HC and Coral. 



number M at selected times, as indicated. At early times (t = 
1 — 10 Myr) all codes agree very well on both the ionization front 
position and its profile. The only modest differences are found in 
the semi-shielded part of the dense gas (xhi = 0.01 — 1), due to 
variations of the treatment of hard photons, and in the low-density 
gas between the clump and the source, where the neutral fractions 
are affected by the slightly different temperatures found by the dif- 



ferent methods. This is confirming the conclusions reached in Paper 
I that with no (or little) gas motions any differences are due to the 
treatment of the energy equation and the hard photons. The hydro- 
dynamic evolution introduces some differences, particularly in the 
I-front position, but the scatter remains small. 

The temperature profiles generally agree in shape and in the 
position of the flow features, the expanding wind and its leading 
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Figure 40. Test 7 (Photoevaporation of a dense clump): Histograms of the gas temperature at times t = 10 and 50 Myrs. 
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Figure 41. Test 7 (Photoevaporation of a dense clump): Histograms of the flow Mach number at times t = 10, and 50 Myrs. 



shock. The main differences are in the amplitude, which varies by 
up to 50%, except for the cold, shielded gas at the back of the re- 
maining dense clump at t = 10 Myr (at position x/Lbox 0.8), 
where the variation between results reached an order of magnitude. 
This large variation does not affect the later-time evolution con- 
siderably, however. The pressure and Mach number profiles (Fig- 
ure 44) show similar trends, with very small differences during the 
early evolution, growing to somewhat larger ones at later times, but 
with all prominent flow features agreeing in both nature and posi- 
tion. 



4.4 Summary and Conclusions 

In this work we compared the results from 10 directly coupled hy- 
drodynamics and radiative transfer codes on three test problems of 
astrophysical interest - H II region expansion in initially uniform 
gas, as well as internal and external photoevaporation of dense 
clumps of galactic-like size and density. Our aims are to validate 
our codes and test their reliability. Our test problems, while cho- 
sen to be relatively simple and clean, nevertheless cover a wide 
range of regimes applicable to photoionization-driven astrophysi- 
cal flows, including propagation of fast (R-type) and slow (D-type) 
I-fronts, shock creation and supersonic photoevaporative winds. All 



the data is available on the Radiative Transfer Comparison Project 
wiki-based website, so future code developers can test their codes 
against our results. 

Overall, the agreement is quite good and all codes are gener- 
ally reliable and produce reasonable results. However the results 
also highlighted some important differences between the meth- 
ods. All participating algorithms track fast, R-type I-fronts well, in 
agreement with the results we obtained in Paper I. We note that this 
is not the trivial statement that we simply reproduce our previous 
static density field results, since in this second Comparison Project 
phase there are several codes which are newly developed (RHID, 
LICORICE, Enzo-RT) and therefore did not participate in Paper I, 
and even the ones which were present then have been further de- 
veloped over the intervening period and are thus not identical to 
the versions used in Paper I. 

Again, as we found in Paper I, the treatment of multi- 
frequency radiative transfer and particularly of the hard tail of the 
photon spectrum varies significantly among the methods and yields 
correspondingly large range of temperature and ionization struc- 
ture just beyond the I-front itself. We showed with a specific ex- 
ample that the spectral energy distribution of the ionizing source 
changes the I-front structure and shocked flow features consider- 
ably. Monochromatic light yields much sharper I-fronts and shocks 
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Figure 42. Test 7 (Photoevaporation of a dense clump): Line cuts of the neutral fraction along the axis of symmetry through the centre of the clump at times 
t = 1 Myr, 10 Myr and 50 Myr (left to right). 




Figure 43. Test 7 (Photoevaporation of a dense clump): Line cuts of the temperature along the axis of symmetry through the centre of the clump at times t = 1 
Myr, 10 Myr and 50 Myr (left to right). 



and certain flow features like the double-peaked profile found in 
Test 5 disappear altogether. 

For static density distributions the variations in the multi- 
frequency radiative transfer treatment had little effect on the I- 
front positions and propagation speeds since those are largely deter- 
mined (apart from recombinations-related effects) by simple pho- 
ton counting and balancing this number against the number of 
atoms entering the front. For dynamically-coupled evolution, this 
changes and there are significant feedback effects, with the radia- 
tive transfer effects affecting the gas dynamics and vice versa. For 
example, pre-heating by hard photons, or lack of it, can affect the 
dynamics significantly. More specifically, higher pre-heating re- 
sults in shocks, e.g. ones typically leading D-type I-fronts, which 
are weaker and faster-propagating, and vice-versa. The internal 
structure of such a front and the relative spacing between the shock 
and the I-front can also change considerably. Shocks created by 
photoheating effects tend to be relatively weak, with Mach numbers 
of a few or less. The density compression resulting from them is 
strongly dependent on the pre-heating by hard photons, but gener- 
ally did not exceed factors of L5-2. The profiles of the fluid quanti- 
ties in supersonically expanding regions (e.g. the photoevaporative 
wind in Test 7) show good agreement among the different methods. 

Significant differences were noted in the numerical diffu- 
sivity of the methods. Numerical diffusion could be due to ei- 



ther the radiative transfer method employed (e.g. the moment 
method OTVET used in HART), or the hydrodynamics (SPH in 
LICORICE). Higher diffusion could have notable effects on some 
properties of the flow (features become smoother, high contrasts are 
diminished), but seems to have modest effects on the overall gross 
features and the basic dynamics remains largely unaffected. How- 
ever, care should be taken when using such methods for problems 
in which the sharp features might matter, e.g. enhanced molecule 
formation due to shocks. 

The propagation of an accelerating I-front down a steep (1/r^) 
density profile proved to be a quite difficult problem and several 
codes developed significant instabilities, while the rest did not. 
While there are a number of physical instabilities which can de- 
velop in similar situations, as we discussed in some detail, in this 
particular case the instabilities we observed proved to be numer- 
ical in nature. The most severe one was the carbuncle instability 
or odd-even decoupling, which in some cases affects low-diffusion 
hydrodynamic solvers (here a Roe Riemann solver). This problem 
can be eliminated by either adding some artificial diffusion or using 
a more diffusive hydrodynamic solver. 

In summary, we have found a considerable level of agreement 
between the wide variety of radiative transfer and hydrodynam- 
ics coupled methods participating in this project. The basic flow 
features and their evolution are reproduced well by all the meth- 
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Figure 44. Test 7 (Photoevaporation of a dense clump): Line cuts of the pressure at times t = 10 Myr (left), and 50 Myr (centre) and of the Mach number at 
time t = 50 Myr (right) along the axis of symmetry through the centre of the clump. 



ods. There are some variations whose origins we did our best to 
understand. The recurring differences were mostly due to the dif- 
ferent treatment of the energy equation and the transfer of multi- 
frequency radiation. There were also some problems specific to 
certain methods which we discussed in detail. While none of the 
codes gave any obviously unphysical or incorrect results and all 
largely agreed with each other, some of the methods were clearly 
less suited for certain problems. No method is universally applica- 
ble to all astrophysical situations and every one of the participating 
codes showed some behaviour discrepant with the majority in one 
respect or another. Care should therefore be taken in applying any 
given algorithm to a new type of problem and detailed testing is 
always advised. 
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